[numcases numdims numbatches]=size(batchdata);

if restart ==1
    restart=0;
    epoch=1;
    
    % Initializing symmetric weights and biases. 
    vishid     = 0.01*randn(numdims, numhid);
    hidbiases  = zeros(1,numhid);
    visbiases  = zeros(1,numdims);

    poshidprobs = zeros(numcases,numhid);
    neghidprobs = zeros(numcases,numhid);
    posprods    = zeros(numdims,numhid);
    negprods    = zeros(numdims,numhid);
    vishidinc  = zeros(numdims,numhid);
    hidbiasinc = zeros(1,numhid);
    visbiasinc = zeros(1,numdims);
end

tmp_sum = 0.0;
count_movie_rate_times = sum(rate_matrix);
for m=1:size(count_movie_rate_times,2)/5
    tmp_sum = sum(count_movie_rate_times(1,m*5-5+1:m*5));
    for n=1:5
        if count_movie_rate_times(1,m*5-5+n)~=0
            visbiases(1,m*5-5+n) = log(count_movie_rate_times(1,m*5-5+n)/(0.0+tmp_sum));
        else
            visbiases(1,m*5-5+n) = 0.0;
        end
    end
end

tic
for epoch = epoch:maxepoch
    rand_index = randperm(numbatches);
    
    for batch = 1:numbatches
        %fprintf(1, 'epoch %4i batch %4i\n', epoch, batch);
        %%%%%%%%% START POSITIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        data = batchdata(:,:,rand_index(1,batch));
        %flag_rate_template = batchdata_real(:,:,rand_index(1,batch));
        
        poshidprobs = 1./(1 + exp(-data*vishid - repmat(hidbiases,numcases,1)));    
        
        posprods    = data' * poshidprobs;
        poshidact   = sum(poshidprobs);
        posvisact = sum(data);
    
        %%%%%%%%% END OF POSITIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        poshidstates = poshidprobs > rand(numcases,numhid);

        %%%%%%%%% START NEGATIVE PHASE%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        negdata = poshidstates*vishid' + repmat(visbiases,numcases,1);
        %negdata = negdata.*flag_rate_template;
        
 %{
        tmpsum = 0;
        for m = 0:num_items-1
            for n = 1:numcases
                tmpsum = max(negdata(n,m*5+1:m*5+5));
                negdata(n,m*5+1:m*5+5) = exp(negdata(n,m*5+1:m*5+5)-tmpsum);
                tmpsum = sum(negdata(n,m*5+1:m*5+5));
                negdata(n,m*5+1:m*5+5) = negdata(n,m*5+1:m*5+5)/tmpsum;
            end
        end
        %}
        for m=1:numcases
            reshap = reshape(negdata(m,:),5, num_items);
            reshap = exp(reshap - repmat(max(reshap), 5, 1));
            reshap = reshap./repmat(sum(reshap), 5, 1);
            negdata(m,:) = reshape(reshap, 1, 5*num_items);
        end
        
        negdata = negdata.*batchdata_real(:,:,rand_index(1,batch));
        
        neghidprobs = 1./(1 + exp(-negdata*vishid - repmat(hidbiases,numcases,1)));  
        
        negprods  = negdata'*neghidprobs;
        neghidact = sum(neghidprobs);
        negvisact = sum(negdata); 

        %%%%%%%%% END OF NEGATIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       
        if epoch>5,
            momentum=finalmomentum;
        else
            momentum=initialmomentum;
        end;

        %%%%%%%%% UPDATE WEIGHTS AND BIASES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
        if batch == numbatches
            count_last = sum(data')>0;
            vishidinc = momentum*vishidinc + ...
                epsilonw*( (posprods-negprods)/sum(count_last) - weightcost*vishid);
            visbiasinc = momentum*visbiasinc + (epsilonvb/sum(count_last))*(posvisact-negvisact);
            hidbiasinc = momentum*hidbiasinc + (epsilonhb/sum(count_last))*(poshidact-neghidact);
        else
            vishidinc = momentum*vishidinc + ...
                epsilonw*( (posprods-negprods)/numcases - weightcost*vishid);
            visbiasinc = momentum*visbiasinc + (epsilonvb/numcases)*(posvisact-negvisact);
            hidbiasinc = momentum*hidbiasinc + (epsilonhb/numcases)*(poshidact-neghidact);
        end
        vishid = vishid + vishidinc;
        visbiases = visbiases + visbiasinc;
        hidbiases = hidbiases + hidbiasinc;

        %%%%%%%%%%%%%%%% END OF UPDATES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    end
   
   % if mod(epoch,10)==0
        MAE
   % end
        toc
        tic

    %fprintf(1, 'epoch %4i error %6.1f  \n', epoch, errsum);
    %error_plot(1,epoch) = errsum;
end;

plot(error_plot);
